setwd( "/data1/bsi/BORA_processing/devel/data_HTHGU/illumina1/MayoGAP/mayogap_emerge/raw/clean/beagle/run_beagle/chr22/1000genome/")
map<-read.table("map.txt",sep="\t",head=T)
segdat<-read.table("segdat.txt",sep="\t",head=T)
outdir <- "/data1/bsi/BORA_processing/devel/data_HTHGU/illumina1/MayoGAP/mayogap_emerge/raw/clean/beagle/run_beagle/chr22/result_new_paper_graph_11_01_2011"
png(file="chr22_median_dosage_diff_1000genome_yrange0.10.png",res=100, width=1600, height=1600)
yrange<-c(0,0.10)
plot(lowess(map$bp, map$meddiff.seed1,f=0.00001), type='l', ylim=yrange, main="Median of Absolute Value of Dosage Difference", xlab="BP", ylab="Median Absolute Difference",lty=2)
lines(lowess(map$bp, map$meddiff.halfmb,f=0.00001),lty=4, col="green")
lines(lowess(map$bp, map$meddiff.seed2,f=0.00001),lty=3, col="red")
lines(lowess(map$bp, map$meddiff.quatermb,f=0.00001), lty=5, col="cyan")
lines(lowess(map$bp, map$meddiff.halfquatermb,f=0.00001), lty=6, col="purple")
lines(lowess(map$bp, map$meddiff.quaterquatermb,f=0.00001), lty=7, col="magenta")
lines(lowess(map$bp[!is.na(map$r2)], yrange[2]-yrange[2]*map$r2[!is.na(map$r2)],f=0.00001), lty=1, col="orange")
abline(v=segdat$seg.halfmb, col="green", lty=4)
abline(v=segdat$seg.quatermb, col="cyan", lty=5)
abline(v=segdat$seg.halfquatermb, col="purple", lty=6)
abline(v=segdat$seg.quaterquatermb, col="magenta", lty=7)
legend("top", col=c("orange","black", "red", "green","cyan","purple","magenta"), lty=c(1,2,3,4,5,6,7), c("Scaled 1-r2,\nPeaks are\nareas of low quality","Different\nSeed1", "Different\nSeed2","0.5MB", "0.25MB"," 0.125MB"," 62.5KB"),horiz=T)
dev.off()
